% The code is based on the following paper:
%   Chen, T. Y., Lin, Y. L., and Tzeng, L. Y., forthcoming.
%   Estimating probability weighting functions through option pricing bounds. Review of Asset Pricing Studies.
%
% Copyright: Tzu-Ying Chen, Yo-Lan Lin, Larry Y. Tzeng
% Date: January 30, 2024

clear

%% Setting 
AllFP = 41;

%% Load Data
FileName = ['AllTradingDate19962021_Monthly.mat'];
load(FileName, 'Target_AllDate');             
clear FileName

FileName = ['StockReturn_19262021_SP500.txt'];
Data_RET = load(FileName);    
Data_RET = sortrows(Data_RET, 1);                                          
clear FileName

%% Estimate the Return Distribution with EGARCH(1,1) Model
HISWindow_EST = 3800;                                                     
HISWindow_EI = 3800;
NumPaths = 250000;

RET_Monthly = nan(size(Target_AllDate, 1), NumPaths);                      
for d = 1:size(Target_AllDate, 1)
    Target_Date = Target_AllDate(d, 1);
    Target_TTM = Target_AllDate(d, 3);

    % Load Data
    FileName = ['OP' num2str(Target_Date) '_TTM' num2str(Target_TTM) '.mat'];
    load(FileName, ...
         'Data', 'Index_Date', 'Index_TTM');
    clear FileName 
    
    Target_Date = Data(1, Index_Date);                                    
    Target_TTM = Data(1, Index_TTM);                                       
    clear Index_Date Index_TTM Data

    % EGARCH(1,1) Model    
    HISWindow = HISWindow_EST;
    Index_HISDate = find(Data_RET(:, 1) < Target_Date);
    Index_HISDate = Index_HISDate((end - HISWindow + 1):end);              
    RET = log(1 + Data_RET(Index_HISDate, end));
    Model_GARCH = egarch('GARCHLags', 1, 'ARCHLags', 1, 'LeverageLags', 1, 'Offset', nan);
    EstModel_GARCH = estimate(Model_GARCH, RET, 'Display', 'Off');
    EstMean = EstModel_GARCH.Offset;                                       
    EstOmega = EstModel_GARCH.Constant;                                    
    EstAlpha = cell2mat(EstModel_GARCH.ARCH);                              
    EstBeta = cell2mat(EstModel_GARCH.GARCH);                              
    EstGamma = cell2mat(EstModel_GARCH.Leverage);                          
    V_GARCH = infer(EstModel_GARCH, RET);          
    I_GARCH = (RET - EstMean) ./ sqrt(V_GARCH);                                              
    clear Index_HISDate HISWindow RET Model_GARCH
 
    HISWindow = HISWindow_EI;
    NumPeriods = max(AllFP);
    Data_EI = I_GARCH((end - HISWindow + 1):end);
    Index_EI = randi(HISWindow, NumPaths, NumPeriods);                     
    I_GARCH_Forecast = Data_EI(Index_EI);
    V_GARCH_Forecast = nan(size(I_GARCH_Forecast));                        
    for i = 1:NumPeriods
        if i > 1
            V_GARCH_Forecast(:, i) = exp(EstOmega + ...
                                         EstAlpha * (abs(I_GARCH_Forecast(:, i - 1)) - sqrt(2 / pi)) + ...
                                         EstBeta * log(V_GARCH_Forecast(:, i - 1)) + ...
                                         EstGamma * I_GARCH_Forecast(:, i - 1));
        else
            V_GARCH_Forecast(:, i) = exp(EstOmega + ...
                                         EstAlpha * (abs(I_GARCH(end)) - sqrt(2 / pi)) + ...
                                         EstBeta * log(V_GARCH(end)) + ...
                                         EstGamma * I_GARCH(end));     
        end
    end
    clear HISWindow 

    RET_GARCH_Forecast = EstMean + sqrt(V_GARCH_Forecast) .* I_GARCH_Forecast;

    % Output
    FileName = ['DG_GrossRET_Monthly_' num2str(Target_Date)];
    for j = 1:length(AllFP)
        FP = AllFP(j);
        RET_Monthly = sum(RET_GARCH_Forecast(:, 1:FP), 2);
        RET_Monthly = exp(RET_Monthly);                                    

        save(FileName, ...
             'RET_Monthly', ...
             'EstModel_GARCH', 'V_GARCH', 'I_GARCH', 'V_GARCH_Forecast', 'I_GARCH_Forecast', 'RET_GARCH_Forecast', ...
             'HISWindow_EST', 'HISWindow_EI', 'NumPaths');         
        clear FP RET_Monthly       
    end
    clear FileName
    
    clear Target_Date Target_TTM NumPeriods
    clear EstModel_GARCH EstMean EstOmega EstAlpha EstBeta EstGamma V_GARCH I_GARCH
    clear Index_EI Data_EI V_GARCH_Forecast I_GARCH_Forecast RET_GARCH_Forecast
    clear RET_Monthly
end
clear Data_RET
clear d i j